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Doping evolution of spin and charge excitations in the Hubbard model 
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To shed light on how electronic correlations vary across the phase diagram of the cuprate su¬ 
perconductors, we examine the doping evolution of spin and charge excitations in the single-band 
Hubbard model using determinant quantum Monte Carlo (DQMC). In the single-particle response, 
we observe that the effects of correlations weaken rapidly with doping, such that one may expect the 
random phase approximation (RPA) to provide an adequate description of the two-particle response. 

In contrast, when compared to RPA, we find that significant residual correlations in the two-particle 
excitations persist up to 40% hole and 15% electron doping (the range of dopings achieved in the 
cuprates). These fundamental differences between the doping evolution of single- and multi-particle 
renormalizations show that conclusions drawn from single-particle processes cannot necessarily be 
applied to multi-particle excitations. Eventually, the system smoothly transitions via a momentum- 
dependent crossover into a weakly correlated metallic state where the spin and charge excitation 
spectra exhibit similar behavior and where RPA provides an adequate description. 


I. INTRODUCTION 

A full understanding of the cuprate phase diagram 
has been prevented in part by difficulties in obtaining 
well-controlled theories for the doping evolution of elec¬ 
tronic excitations. This is due to the lack of exact solu¬ 
tions to the Hubbard model in two dimensions, the stan¬ 
dard model thought to contain the low-energy physics of 
the cuprates that captures aspects of magnetic proper¬ 
ties seen in experiments^. At low doping, strong cou¬ 
pling treatments of the large-U Hubbard model describe 
the dispersion and intensity of magnon excitations ob¬ 
served via neutron scattering^. At high doping, it is be¬ 
lieved that spin excitations behave like weakly interacting 
particle-hole excitations governed by the underlying free 
particle kinetic energy, with a weak influence from the 
Hubbard U. If correct, this limit can be adequately rep¬ 
resented by the random phase approximation (RPA) as a 
proxy for more exact treatments^. Many studies have as¬ 
sumed weak correlations in doped cuprates, so RPA has 
been used to address spin and charge excitations observed 
by neutron and Raman scattering, as well as the forma¬ 
tion of a d-wdive superconducting ground state^’^ 

However, a set of surprising results have emerged from 
recent resonant inelastic x-ray scattering measurements 
(RIXS) in a variety of cuprates^^. In the hole-doped 
La2— xSr3^Cu04, YBa2Cu306-|-3^, 

Tl2Ba2Cu06+cc, and Bi2Sr2CaCu208+a: families^^^^^, 
high-energy magnons or paramagnons on the antiferro¬ 
magnetic zone boundary (AEZB) persist from the parent 
compounds into the heavily overdoped regime, showing 
little doping dependence up to 40 % hole doping where 


the system is believed to exhibit Eermi-liquid-like behav¬ 
ior in the single-particle response or transport. In the 
electron-doped Nd2-ccCea,Cu04 family^^’^^, momentum- 
dependent low-energy charge excitations have been found 
over a large energy range and surprisingly, the magnetic 
excitations harden to high energies when doped beyond 
the antiferromagnetic phase. Both computational and 
analytical techniques have been brought to bear on these 
results. Exact diagonalization and determinant quantum 
Monte Carlo (DQMC) have captured the momentum and 
doping dependence of the AEZB paramagnons^^. On the 
other hand, RPA has been used to conclude that the 
collective spin and charge excitations have similar low- 
energy behavior^^. Beyond the spin response, the charge 
excitations may be associated with charge ordering ob¬ 
served in the cuprate pseudogap regime, so it is useful to 
systematically explore how both spin and charge excita¬ 
tions evolve with doping throughout the Brillouin zone, 
and to what degree the response functions can be ap¬ 
proximated using RPA. 

In this study, we compute spin and charge susceptibil¬ 
ities of the single-band Hubbard model^’^^’^^ throughout 
the first Brillouin zone (Eig. I) via DQMC, a numerically 
exact imaginary-time auxiliary-field technique^^~^^. The 
susceptibilities are compared to those calculated with 
rPA 26 ^ ^ formalism originally developed for weakly inter¬ 
acting systems that is expected to become an increasingly 
good approximation as the doping level increases. Our 
calculations reveal that the influence of correlations on 
the multi-particle excitations persists to higher dopings 
than suggested by quantities related to the single-particle 
response^^. At even higher doping, RPA provides an ad- 
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FIG. 1: Cartoon of the first Brillouin zone showing the an¬ 
tiferromagnetic zone boundary (AFZB) and high-symmetry 
cuts. The nodal direction, from zone center to (tt, tt), is indi¬ 
cated with a red arrow. The antinodal direction, from zone 
center to (tt, 0), is marked with a green arrow. 


equate description of the response functions, which show 
a smooth momentum-dependent crossover into a weakly 
correlated metallic state. 


II. MODEL AND METHODS 

The single-band Hubbard Hamiltonian describes 
strongly correlated electrons on a lattice: 

^ = Y. + h-C.) + Y t'iMaCja + h-C.) 

{hj)cr {i,j}cr 

-f^Y^i'^ + uY ni'friii ( 1 ) 

ia i 


where cj^ (Ca) creates (annihilates) a particle with spin 
a on site i, and rii^ = is the number opera¬ 

tor. The nearest- and next-nearest-neighbor hoppings are 
controlled by t and t', respectively, and U is the on-site 
Coulomb interaction that penalizes double occupancy. 
We work with a canonical parameter set, t' = — 0.3t and 
U = 8t, so the ground state is a strongly correlated Mott 
insulator in the undoped system^^ and upon hole-doping 
the system will possess a hole-like Fermi surface. As 
usual all energies are expressed in units of t. The chemi¬ 
cal potential /i is adjusted to give between 15% electron 
and 75% hole doping. Without particle-hole symmetry, 
DQMC exhibits a significant sign problem^^, so we work 
at an inverse temperature /3 = 3/t to give reasonable 
statistics. The imaginary-time spin and charge correla¬ 
tors are computed as 


Xs,c(q, t) = r)dtc(q, 0)), (2) 

where Os = Oc = Sj e*'*'*^’(nq + 

n^), and riia- is the number operator. The correlators 
are analytically continued using the Maximum Entropy 
method to obtain the real-frequency susceptibilities^^. 

RPA susceptibilities are determined for the same tem¬ 
perature and doping range as DQMC, and are normalized 





FIG. 2: The spin susceptibilities along high-symmetry cuts 
in the Brillouin zone are calculated using DQMG [panels 
(al)-(a3)] and RPA [panels (bl)-(b3)], for three hole dopings 
[n=1.00 in (al) and (bl), n=0.85 in (a2) and (b2), and n=0.60 
in (a3) and (b3)]. The maximum of the color scale is set by the 
highest intensity in each panel, (c) The peak energies in (al)- 
(a3) and (bl)-(b3) are plotted versus momentum to highlight 
the doping trends. Spin susceptibilities at the representative 
momentum points (7r,7r) and (tt, 0) from DQMG [(dl) and 
(d2), respectively] and RPA [(d3) and (d4), respectively] are 
shown for hole doping ranging from 0% to 40%. 


to the DQMC two-particle sum rule for ease of compari¬ 
son. The RPA susceptibilities are computed from: 


X 


RPA 

s,c 


(q,w) 


xo(q,^) 

1 =F exo(q,w)’ 


(3) 


where — (+) corresponds to the spin (charge) susceptibil¬ 
ity, U is the effective interaction strength, and xo is the 
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non-interacting Lindhard susceptibility^’^^ 


xo(q,cj) 


/(^fe+g) — fi^k) 
^ k ~ - e/c) 


( 4 ) 


The bandstructure is Ck = —2t{coskx + cosky) — 
At' cos kx cos ky and the Fermi function is f{ek) = {1 + 
exp [/d(e/c — with the chemical potential /i control¬ 

ling the filling. The inset in Fig. 2(d3) shows Umax = 
Rex^^ 0), the location of the new pole in plotted 

versus filling for (tt, tt), as well as for (tt, 0) and (7r/2,7r/2) 
on the AFZB. U is most tightly constrained by (7r,7r). 

To isolate doping-dependent bandstructure effects 
from doping-dependent changes to effective interactions 
that can be harder to predict, we use a single value of 
U to calculate the RPA susceptibility for all momenta 
and all doping levels. Because the low-doping behav¬ 
ior of the spin response in cuprates, and the single-band 
Hubbard model, is well known, we set U = 3t such that 
the RPA response mimics the spin response near the AF 
instability at (tt, tt) for 0% doping which captures the ex¬ 
perimental neutron scattering results^^ and matches the 
DQMC trends by construction. These low-energy spin 
excitations near (tt, tt) also may be integral to pairing in 
a broad class of superconductors^, and their doping evo¬ 
lution will thus be well captured by both computational 
methods. We note that this value for U falls within the 
range of values obtained in previous studies^’^^O- how¬ 
ever, one should note again that for the chosen set of pa¬ 
rameters the undoped, single-band Hubbard model pos¬ 
sesses a Mott insulating ground state, and should not be 
viewed as weakly correlated. The expression for x^^^ 
does not contain a new pole, so its behavior is primar¬ 
ily determined by that of xo and thus set by the non¬ 
interacting bandstructure. 


III. DOPING EVOLUTION OF SPIN AND 
CHARGE SUSCEPTIBILITIES 

Figure 2 summarizes the hole doping evolution of the 
DQMC [(al)-(a3)] and RPA [(bl)-(b3)] spin susceptibil¬ 
ities along the nodal (from q = (0,0) to q = (7r,7r)) and 
antinodal (from q = (0, 0) to q = (tt, 0 )) cuts in the first 
Brillouin zone. For momenta near (tt, tt), doping both re¬ 
duces and shifts spectral weight to higher energies, with 
a good comparison between DQMC and RPA as doping 
increases, a mirror of the results from neutron scatter¬ 
ing experiments^^and a by-product of the choice of 
U. On the other hand, along the antinodal direction the 
peak in the spin response remains unchanged with dop¬ 
ing in the DQMC calculations but softens considerably in 
RPA, highlighted in Fig. 2(c). The strong disagreement 
between DQMC and RPA in large portions of the Bril¬ 
louin zone as doping increases indicates that the infiuence 
of correlations on the spin response remains considerable 
to relatively large doping levels. 

Since (tt, tt) and (tt, 0 ) are representative of the behav¬ 
ior along the nodal and antinodal directions, respectively. 





FIG. 3: The charge susceptibilities along high-symmetry cuts 
in the Brillouin zone are calculated using DQMC [panels 
(al)-(a3)] and RPA [panels (bl)-(b3)], for three hole dopings 
[n=1.00 in (al) and (bl), n=0.85 in (a2) and (b2), and n=0.60 
in (a3) and (b3)]. The maximum of the color scale is set by the 
highest intensity in each panel, (c) The peak energies in (al)- 
(a3) and (bl)-(b3) are plotted versus momentum to highlight 
the doping trends. Charge susceptibilities at the represen¬ 
tative momentum points (tt, tt) and (tt, 0) from DQMC [(dl) 
and (d2), respectively] and RPA [(d3) and (d4), respectively] 
are shown for hole doping ranging from 0% to 40%. 


Fig. 2(dl)-(d4) focuses on the evolution from 0% to 40% 
hole doping at those momenta. At 0% doping, the (tt, tt) 
spin response peaks strongly at low energy, but it broad¬ 
ens and hardens with increasing doping, as determined 
by the non-interacting bandwidth. These findings agree 
with previous DQMC and dynamical cluster approxima¬ 
tion calculations^^”^^. Due to proximity to the pole in 
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FIG. 4: The spin and charge susceptibilities are calculated by 
DQMC and RPA for different levels of electron doping. While 
the DQMC spin susceptibility shows the same doping trend as 
RPA at (7r,7r) and (7r,0), suggesting decreasing correlations, 
the RPA charge susceptibility is dominated by the band edge, 
unlike DQMC. 


FIG. 5: The spin and charge susceptibilities calculated by 
DQMC and RPA are shown for representative momenta at 
55% and 75% hole doping to illustrate the gradual weakening 
of correlations at extremely high doping. 


Eq. 3, the low-doping RPA spin excitations at (tt, tt) gen¬ 
erally occur at lower energies than in DQMC and their 
peaks are much sharper. In addition, the RPA response 
has longer tails set by the full non-interacting bandwidth. 
At (7r,0), the DQMC and RPA spin responses remain 
quite different up to at least 40% hole doping. The 
RPA spin susceptibilities soften significantly with dop¬ 
ing, whereas the DQMC spin response shows persistent 
excitations in agreement with experiments. 

Figure 3 summarizes the hole doping evolution of the 
charge susceptibilities as calculated using DQMC [(al)- 
(a3)] and RPA [(bl)-(b3)] throughout the first Brillouin 
zone. Along both the nodal and antinodal cuts, the spec¬ 
tral weight of the response evaluated using DQMC is lo¬ 
cated at high energies (determined at 0% doping by the 
Hubbard U) and the spectra show a charge gap that de¬ 
creases with doping^^’^^. On the other hand, the RPA 
susceptibility shows no charge gap and is dominated by 
the peak in the Lindhard susceptibility at the band edge. 
With doping, the lineshape of the RPA charge response 
along the nodal cut changes little [Figs. 3(bl)-(b3)], al¬ 
though spectral weight along the antinodal direction in¬ 
creases at lower energies in agreement with the response 
from DQMC. Figure 3(c) highlights the doping evolution 
of the peak positions and Figs. 3(dl)-(d4) focus on a 
more detailed doping dependence for representative mo¬ 
menta in the nodal, (7r,7r), and antinodal, (7r,0), direc¬ 
tions. 

Upon electron doping, the DQMC and RPA spin sus¬ 
ceptibilities (Fig. 4) show similar behavior, hardening 
and broadening with increasing doping at (7r,7r), again 
by construction. At (7r,0), both Imx^^^ and Imx^^^^ 
harden and decrease in intensity with doping. The spin 
response thus exhibits an electron-hole doping asymme¬ 
try, as the doping trend near (tt, 0) is different from 


that on the hole-doped side (Fig. 2). This suggests 
that correlation effects on the spin susceptibilities be¬ 
come less relevant more rapidly with electron than hole 
doping. In the charge channel, the DQMC and RPA 
responses exhibit different doping trends similar to the 
evolution on the hole-doped side: the charge gap along 
the nodal and antinodal directions at 0% doping in the 
DQMC charge susceptibility vanishes with electron dop¬ 
ing whereas the RPA charge susceptibility shows no gap, 
and the low-energy structure in Imx^^^ remains essen¬ 
tially unchanged up to 15% electron doping for both 
(tt, tt) and (tt, 0), again manifesting the system’s electron- 
hole doping asymmetry. 

As the system is doped to extremely high levels, 
the correlations continue to weaken until the system 
smoothly transitions into a weakly correlated metal, 
where RPA provides a fairly adequate description of the 
spin and charge responses across the Brillouin zone. As 
shown in Fig. 5, at 55% and even more so at 75% 
hole doping, Imx^^^^ and Imx^^^ agree both quali¬ 
tatively, and even to some degree quantitatively, at low 
energies. In addition, the DQMC spin peak broadens sig¬ 
nificantly with increasing doping until it closely resembles 
the charge susceptibility at 75% hole doping, indicating 
that the response in the spin channel essentially can be 
viewed as a trivial spin flip on top of the charge excita¬ 
tions. Figure 5 thus demonstrates that the influence of 
correlations on the multi-particle response can decrease 
significantly, but only at doping levels well beyond those 
where weak coupling approaches typically already have 
been applied for the cuprates. 
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IV. CORRELATIONS IN SINGLE- VERSUS 
TWO-PARTICLE CORRELATORS 


An important subtlety in examining how correlations 
evolve with doping is that single- and multi-particle 
quantities exhibit fundamentally different renormaliza¬ 
tions. The spin and charge susceptibilities reveal that 
correlations in two-particle quantities persist to signifi¬ 
cantly higher dopings than suggested by single-particle 
correlators^^, necessitating the study of both to resur¬ 
rect the full richness of the phase diagram. Figure 6(a) 
shows the doping dependence of the compressibility for 
t' = 0 and —0.3t, with the next-nearest-neighbor hopping 
breaking particle-hole symmetry. Because the compress¬ 
ibility is zero in the Mott plateau and finite elsewhere, it 
could be interpreted as a proxy of how strong the corre¬ 
lations are. For both = 0 and —0.3t, the compressibil¬ 
ity vanishes near half filling, corresponding to the Mott 
plateau, as expected. However, it recovers rapidly with 
both hole and electron doping, suggesting that the cor¬ 
relation effects weaken rapidly away from half filling. 

On the other hand, the nearest-neighbor equal-time 
spin-spin correlation function [Fig. 6(b)] suggests that 
correlations extend to higher doping levels. For = 0, 
the spin correlations persist well away from half fill¬ 
ing with both hole and electron doping. Next-nearest- 
neighbor hopping strongly suppresses the correlation 
function above 30% hole doping, suggesting that corre¬ 
lations have weakened significantly by that doping level, 
but the magnitude of the correlation function is actually 
enhanced slightly on the electron-doped side. The equal¬ 
time spin-spin correlation function thus implies that cor¬ 
relation effects persist to significantly higher doping levels 
than can be seen in the compressibility. 

The difference between the double occupancy D = 
from its uncorrelated value of ((n)/2)^ provides a 
third way of studying how far correlations extend away 
from half filling in the charge response [Fig. 6(c)]. Due to 
its local nature, this quantity is more sensitive to corre¬ 
lations. Like the spin-spin correlation function, it shows 
that they weaken more slowly with doping than the com¬ 
pressibility suggests. However, unlike the spin-spin cor¬ 
relation function, ((n)/2)^ — D shows little particle-hole 
asymmetry in the doping trend even with next-nearest- 
neighbor hopping, implying a fundamental difference be¬ 
tween the spin and charge responses. 

Figure 6 also highlights the particle-hole symmetry¬ 
breaking effect of next-nearest-neighbor hopping on 
single- and multi-particle quantities. Although the com¬ 
pressibility shows little doping dependence above ^ 15% 
hole or electron doping when t' = 0, it exhibits notice¬ 
able doping dependence on the hole-doped side when 
t' = —0.3t. The spin correlations are suppressed on the 
hole-doped side by next-nearest-neighbor hopping, while 
they are enhanced on the electron-doped side. Although 
D is barely impacted by t' = — 0.3t because of the dou¬ 
ble occupancy’s local nature, single- and multi-particle 
quantities that are sensitive to longer-range hopping will 
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FIG. 6: The doping trends of the (a) compressibility, (b) 
equal-time spin-spin correlation function, and (c) difference 
between the double occupancy D = horn its uncor¬ 

related value of ((n)/2)^ are shown for next-nearest-neighbor 
hopping F = 0 and —0.3t. 


display significant particle-hole doping asymmetry. 


V. CONCLUSIONS 

In this study, we have investigated spin and charge sus¬ 
ceptibilities in the single-band Hubbard model to under¬ 
stand the influence of correlations on the multi-particle 
response as a Mott insulator evolves into a weakly cor¬ 
related metal. The naive expectation based on probes 
more sensitive to single-particle properties has been for 
correlations to weaken rapidly with doping such that 
the cuprates cross-over to a more Fermi-liquid-like be¬ 
havior for overdoping near the edge of the supercon¬ 
ducting dome. The compressibility, a quantity associ- 
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ated with the single-particle response, computed using 
DQMC reflects this behavior as it quickly becomes non¬ 
zero and even saturates away from 0% doping. How¬ 
ever, in addition to the dynamical response functions dis¬ 
cussed here, two-particle equal-time quantities computed 
using DQMC also show significant residual correlations 
to much higher doping levels: the behavior of the spin- 
spin correlation function counters the naive expectation 
that the response in the spin channels represents a simple 
spin flip on top of the charge background, and the differ¬ 
ence in the double occupancy from its uncorrelated value 
provides a direct measure of the residual correlations in 
the system, which extend to relatively high doping lev¬ 
els as concluded from the multi-particle response. When 
compared to RPA calculations, the DQMC-derived spin 
and charge response functions show qualitative differ¬ 
ences that persist across large portions of the Brillouin 
zone and throughout the doping range relevant to the 
cuprates, especially in the charge channel, attributable 
to distinctions between the doping dependence of corre¬ 
lation effects at the single- and multi-particle level. Only 
when doped to extremely high levels will the two-particle 
response represent a system in a weakly correlated metal¬ 
lic state. These conclusions help to elucidate evolution 


away from the Mott insulating ground state and demon¬ 
strate that strong correlations can extend over a larger 
region of the cuprate phase diagram than has been ap¬ 
preciated previously. 
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